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ABSTRACT 

The modeling of damped signals as the impulse response of a pole-zero system 
is considered for a broad range of pole-zero modeling algorithms. The goal is to 
obtain the best possible fit between the model impulse response and the modeled sig- 
nal. Prony’s method, the least squares modified Yule- Walker equations (LSMYWE), 
iterative prefiltering, and the Akakie maximum likelihood estimator are compared 
on known test sequences for a variety of model degrading situations (e.g., additive 
noise) to develop an understanding of which methods are most suitable for mod- 
eling real world signals. A correlation domain version of interative prefiltering is 
also introduced. The most robust algorithms are determined to be LSMYWE using 
singular value decomposition and iterative prefiltering (including the correlation do- 
main version). Modeling several laboratory generated short duration acoustic signals 
confirmed the robustness of LSMYWE and iterative prefiltering. It is shown that 
correlation domain iterative prefiltering outperforms standard iterative prefiltering 
when large model orders are required for accurate modeling. Shank’s method was 
determined to be the most effective method of determining the zeros of a pole-zero 
model when a time domain match is required. 
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I. INTRODUCTION 



A. RATIONAL MODELING OF TIME SERIES DATA 

The need to find a compact parametric representation for time series data arises 
in many fields of study. In electrical engineering, finding such representations, or 
models, appears under such topics as optimum control, optimum filtering, system 
identification, model order reduction, waveform encoding, and spectrum estimation. 
Other fields refer to such modeling as time series analysis or forecasting. In digital 
systems, the model parameters take the form of difference equation coefficients or, 
alternatively, transfer function coefficients. The most general form of a difference 
equation is one that uses both feed-forward and feedback coefficients. The corre- 
sponding transfer function is a ratio of two polynomials in the complex variable z 
and is refered to as a rational, or pole-zero, model. 

Historically, the more general pole-zero model has been used in relatively few 
applications compared to all-pole (feedback only) and :ii! zero (feed forward only) 
models. In some cases, an all-pole or all-zero model is llic most appropriate. More 
often, however, the all-pole or all-zero model is chosen because the optimal model 
estimation procedures are better understood and easier to implement than those for 
pole-zero modeling. This is particularly true for the all-pole case which, in many 
situations, can be determined by solving a set of linear equations. Two recent devel- 
opments have led to increased activity in applying pole-zero models: (1) technological 
advances in digital hardware have dramatically reduced computation costs and, (2) 
a greater variety of efficient techniques for estimating pole-zero model parameters is 
now available. 
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Literally hundreds of pole-zero modeling estimation algorithms and applications 
have been published. The majority of these are built around probablistic or stochas- 
tic modeling techniques. This is because stochastic modeling is usually the most 
appropriate to forecasting [Ref. 1, 2, 3] and spectrum estimation [Ref. 4, 5, 6] where 
very little is known about the system input which produced the time series being 
modeled. A deterministic methodology has also been used in which the input and 
output time series are available and the linear system which ‘best’ (usually in a least 
squares sense) produces this cause and effect is determined. This approach is usually 
found under the topics system identification [Ref. 7, 8, 9, 10] and waveform encoding 
[Ref. 11). Another large body of literature which exists in parallel to stochastic and 
deterministic modeling is that of reduced-order modeling [Ref. 12] which largely deals 
with the system control applications of pole-zero modeling. 

Although stochastic and deterministic pole- zero modeling have the same goal 
and use the same mathematical techniques, the distinction is significant in the way 
data is treated and in the performance criterion postulated. Specifically, stationarity 
requirements and assumptions about the probability density functions of random pro- 
cesses in stochastic modeling can be limiting when dealing with real world systems. In 
contrast, deterministic techniques make no specific assumptions about the time series 
to be modeled except, of course, that the form of the model chosen is appropriate. 

One particular deterministic modeling problem that has received little detailed 
attention is that of finding a pole-zero model when the time series being modeled 
is a transient, ‘impulse response’-like waveform. Examples of situations where such 
models may be useful include in modal or shock analysis of mechanical systems [Ref. 
13, 14], antenna response to electromagnetic pulses [Ref. 15], and wavelet estima- 
tion in siesmic studies [Ref. 16]. Pole-zero models make sense for transient wave- 
forms because such a model is structural for many transient signals. That is, impulse 
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response like transients are usually the result of a system of oscillators which have 
been excited for a very short time period relative to the natural frequencies of the 
oscillators. Pole pairs in pole-zero models correspondingly represent the resonant fre- 
quencies of a digital system. System zeros allow the initial conditions or phasing of 
the resonant frequencies to be modeled. This combination of poles and zeros allows 
signals to be modeled based on time domain matching. When an effective time do- 
main match is achieved many concerns about assumptions during modeling become 
moot; an effective time domain match has, by definition, effectively characterized the 
signal in question. Previous work has concentrated on finding only transient model 
poles [Ref. 17, 18, 19, 20] or concentrated on one particular technique of pole-zero 
impulse response matching [Ref. 21]. 

B. THESIS OUTLINE 

This thesis provides a performance comparison of several pole-zero modeling 
procedures applied to the problem of model estimation from impulse response data. 
The procedures compared are chosen to form a cross section of optimality and compu- 
tational complexity of available techniques. In Chapter II, the modeling procedures 
selected for study are described in detail. The remaining three chapters are concerned 
with comparing the performance of these modeling techniques and are organized as 
follows: 

1. To study the specific modeling properties of each method, test impulse response 
sequences are modeled in Chapter Three. Test sequences are constructed to 
simulate the degradations likely to be present in ‘real world’ transient signals 
(e.g. noise, unknown model order). 

2. Using the results obtained in Chapter III, laboratory generated acoustic tran- 
sient data is modeled in Chapter IV. This data is considered to determine the 
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performance of various techniques when modeling the highly complex data char- 
acteristic of real world sources about which very little is usually known. 



3. Chapter V summarizes the main conclusions drawn from the results in Chapters 
III and IV. Recommendations for further study are also presented. 
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II. POLE-ZERO MODELING 



A. OVERVIEW — THE VARIETY OF TECHNIQUES 

Choosing a pole-zero modeling technique from among the many available tech- 
niques can be difficult. The complexity, applicability, and demonstrated effectiveness 
of the different methods are not always well documented. Additional consideration of 
the many refinements that often evolve cis a technique is applied to different problems 
can lead to a perplexing array of tools with which to attack the modeling problem. 
For the particular case of modeling by impulse response matching, little work has ap- 
peared which sorts out the strengths and weaknesses of available modeling techniques 
or, in fact, indicates which methods may or may not be applicable. 

The basic scheme for fitting a pole-zero system impulse response to a given data 
sequence, x(n), is illustrated in Figure 2.1a. This is sometimes referred to as the direct 
model. As we shall see, the formulation of this problem leads to a set of nonlinear 
equations which require the use of iterative techniques to solve. To overcome the 
complexities inherent in solving nonlinear equations, the impulse response matching 
problem can be reformulated as shown in Figure 2.1b. This may be referred to as the 
indirect method. Note that these two formulations are not equivalent: the error of 
the direct method of Figure 2.1a is e<f(n) = x(n) — h(n) while the error for the indirect 
method of Figure 2.1b is e,(n) = 6(n) — a(n)*x(n), where k(n) is the impulse response 
of the pole-zero system being found, b(n) is the corresponding sequence of numerator 
coefficients and a(n) is the corresponding sequence of denominator coefficients. The 
solution of the indirect problem is considered suboptimal in the sense that, except 
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x{n) 





(b) 



Figure 2.1: Two pole-zero modeling problem formulations: (a) direct for- 
mulation and (b) indirect formulation. 
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when the modeling error goes to zero, the effect of minimizing ej(n) is not the same 
as minimizing e<f(n) in the direct method. 

One way to organize pole-zero modeling techniques is shown in Figure 2.2. In 
deterministic waveform matching, the ideal equations relating the proposed model to 
the available input and output data sequences are formed. The solution which ‘best’ 
satisfies these equations is chosen. In this context, ‘best’ usually means the solution 
which minimizes the sum of squares of the equation error. These are essentially the 
Prony type methods [Ref. 17, 22] (when the input is an impulse response) and least 
squares system identification [Ref. 7, 8, 23] methods (for generail input sequences). 
A number of iterative techniques for waveform matching have also been proposed. 
Waveform fitting error [Ref. 24, 25] and inverse filtering error [Ref. 26, 27] are the 
criteria most used. 

Linear stochastic pole-zero modeling techniques rely primarily on estimates of 
second order statistics (auto- and cross-correlations) to estimate model parameters. 
Spectral estimation has been a driving force for these methods which are based on 
solving some form of the modified Yule- Walker equations (see Chapter III). Other 
methods which utilize reflection coefficients [Ref. 28] and higher order statistics [Ref. 
29, 30] have also appeared. 

The maximum likelihood technique seeks parameter estimates for which the ob- 
served data is the most probable in the sense that its conditional probability density 
function (likelihood function) is maximized. This technique is considered to be sta- 
tistically optimum but is quite difficult to use because its implementation generally 
requires the minimization of a highly nonlinear function [Ref. 1, 31, 32, 33]. Four 
widely applied methods from each of the categories of Figure 2.2 will be used in 
this thesis. A number of improvements which have subsequently been suggested for 
these methods will also be considered. Most recent work in pole-zero modeling and 
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spectrum estimation has occured at the internal boundaries of Figure 2.2, i.e. equiv- 
alent linear techniques are sought that perform as well as modeling formulations that 
require solving nonlinear equations [Ref. 34, 35, 28, 36, 37]. The application of 
these newer techniques to transient modeling will not be considered here since we 
expect that the performance of the methods chosen will in most cases bracket the 
performance of these newer techniques with the possible sacrifice of computational 
efficiency. The four pole-zero modeling techniques chosen will be described in the re- 
mainder of this chapter. The transient modeling performance of these methods that 
will then be compared in Chapters III and IV. 



Linear 



Iterative 



Deterministic 



Stochastic 



Equation Error 
Methods 


Waveform Matching 
Inverse Filtering 


Correlation Equation 
Error Methods 


Maximum 
Likelihood Methods 



Figure 2.2; One way of organizing the various pole-zero modeling tech- 
niques. 



B. THESIS MODELING TECHNIQUES 
1. Prony’s Method 

One of the best known indirect techniques for matching a waveform to 
the impulse response of a linear time-invariant system is Prony’s method [Ref. 22]. 
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This method is in fact a special case of least squares system identification in which 
the system input sequence is taken to be a unit impulse and the numerator and 
denominator coefficients are determined separately. 

The pole-zero modeling problem is formulated as follows. The time domain 
difference equation for a general feedback, feed-forward system can be written 

J3a_,x(n -i) = 536iu(n - i) (2.1) 

j=0 i=0 

where n(n) is the input sequence and x{n) is the output sequence. When the input 
sequence is taken to be a unit impulse (unit sample function) and the output is taken 
to be the corresponding impulse response, (2.1) can be expressed in the form of a 
matrix equation, 

■ x(0) 0 

x(l) i(0) 

.x(A^-l) x(N-2) 

where N is the number of data points used and, without loss of generalization, cq is 
set equal to one. 

Equation (2.2) can be solved by partitioning. 



0 

0 

x{N-P) 





■ 1 ■ 




Cl 




Ci . 

1 



Oo 



bQ 

0 



( 2 . 2 ) 



Xb 



a = 



b 

0 



( 2 . 3 ) 



where a = [1 oi • • • ap]^, b = [6o 6i • • • bq]^ and and Xp are the corresponding 
lower and upper partitions of the data matrix in (2.2). The upper partition consists 
of the first Q + 1 rows of the data matrix and the lower partition is composed of the 
remaining rows. 

The solution can then be obtained by first solving the lower partition, 



X^a — 0, 



( 2 . 4 ) 
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for a and then finding b from the upper partition, 

Xsa = b. (2.5) 

l{N = P-\-Q-i-l then (2.2) may have a unique solution and the model impulse 
response will exactly match x(n) forO<n<P + Q4-l. This is referred to as the 
Fade approximation [Ref. 22]. 

In most circumstances, however, the length of the available data sequence 
far exceeds P + Q + 1. It is then desirable to use all available data in setting up (2.2). 
This leads to an overdetermined set of linear equations for the lower partition. No 
exact solution to (2.4) usually exists in this case. The relationship in (2.4) becomes 



X^a = e. 



( 2 . 6 ) 



where e is the equation error that will be present. The solution of (2.4) and (2.6) 
requires the partitioning of as follows, 

X>i = [ ] (2.7) 

where is the first column of X^. If the remaining matrix X^ is square and of full 
rank, then the solution to (2.4) is given by 



a'=(X:,)-^x^ 



( 2 . 8 ) 



where 



1 



a = 



(2.9) 



Otherwise the least squares error solution of (2.6), which minimizes the squared error 
e^e, is given by 

a' = X'+x^ (2.10) 
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where is the Moore-Penrose pseudoinverse of X'^. When X^ is of full (column) 
rank P, the psuedoinverse is given by 



a' = 



( 2 . 11 ) 



(See [Ref. 4, pp. 28-33] for an example of this derivation). Otherwise, the pseudoin- 
verse is defined by 



w T 
^ u/ X>iV, 



( 2 . 12 ) 



(=1 <’■ 

where <r,-, u,-, and v,- are defined by the singular value decomposition (SVD) represen- 



tation of X^ , 



w 



X'x = 2Z<t,u,V( 



(2.13) 



.=0 

The parameters cr,- are the singular values of X^, u,- and v, are the corresponding 
left and right singular vectors, and W is the rank of X^. See [Ref. 38, Ch. 12] for 
a more detailed explanation of singular value decomposition. Once a is known, the 
upper partition of (2.2) can be solved by simply carrying out the matrix multiplication 
described by the left hand side of (2.5). 

2. Modified Yule- Walker Equation Methods 

A well known class of pole-zero modeling techniques is based of solving 
some form of the Modified Yule- Walker Equations (MYWE). These equations can be 
developed by multiplying (2.1) by x{n — k) and taking the expectation of both sides. 
This yields. 



P Q 

Y, ^kTxxin -k) = Y - i) (2.14) 

Jfc=0 t=0 

where rxx{l) is the autocorrelation sequence of the system output and r^xii) is the 
crosscorrelation of the system input and output. If the original input to (2.1) is 
assumed to be a unit variance white noise sequence then the cross correlation, r„i(f), 
is given by 
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(2.15) 



jB[u(n)a;(n — /)] = E\u{n) h{k)u{n — I — k)] 

k 

= Y;^h{k)6{i+k) 

= h{-l). 



Equation (2.14) can be then be written as 



P Q 

0'krxx{n -k) = Yl - ”)• 

fc=o «=o 



or in matrix form, 



^xi(O) 

'I' XX ( 1 ) 



rxx(9). 
^xx(9 + i) 



7’xx( 1) 

^xx(O) 

rxx(9 - 1) 

^xx ( 9 ) 



rxx{-P) 
Txx{-P + 1) 

rxx(9 - P) 

pxx(Q + i — p) 



rxx{Q-\-P) ^xx(Q + P-l) 

~^xx{Q + P + 1) ^xx(Q + P) 



1 

Cl 



1'xx{Q + 1) 

^xx(Q) 

T.%obih{i) 
Ef=i b,h{i - 1) 

6q/i(0) 

0 



(2.16) 



(2.17) 



0 

As in Prony’s method, the solution for the Ck and 6,- coefficients can be 
determined separately. Taking a lower partition of the last P equations in (2.17) 
results in the matrix equation, 

R^a = 0. (2.18) 

where the theoretical values of the elements in are replaced by estimated values. 

If P + 1 autocorrelation lags are used in constructing R^, then (2.18) can 
be solved directly for a. However, if additional reliable lag information is available, we 
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will again desire to extend by letting the index n in (2.16) run beyond Q + P + 1 
resulting in additional equations. This leads to a an overdetermined set of equations, 



R^a = e (2.19) 

which will not in general be satisfied with zero error. As before, application of the 
psuedoinverse results in a least squares solution of (2.19). 

To understand how the MYWE methods can be used to match a time 
series to the impulse response of a pole-zero system observe that (2.16) describes the 
relationship depicted in Figure 2.3. This operation can be equivalently expressed as 

fxx{Ti) = h{n) * h{—n). (2.20) 



If the signal to be modeled is assumed to be a pole-zero system impulse response, then 
for the purpose of implementing (2.17), the signal being modeled can be substituted 
for h{n) in (2.20). 



h{-n) 




r’xx(^) 



) 



Figure 2.3: The system relationship described by (2.16). 



The equations of the upper partition of (2.17) (the first Q equations) are 
not linear in the b coefficients and are generally not solved directly from (2.17). Once 
the denominator coefficients have been determined from (2.18) or (2.19), any of a 
number of techniques are available for finding the desired transfer function numerator 
coefficients. One method already discussed is to set up and solve the upper partition of 
Prony’s method in (2.5). Three other techniques are spectral factorization, Durbin’s 
method, and least squares identification, each of which are described below. 
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a. Spectral Factorization 

Equation (2.16) describes the time domain relationship 

r,,{l)*a{l) = b{l)*h{-l) (2.21) 

where a{l) and b{l) represent the sequences of denominator and numerator transfer 

function coefficients, respectively. Taking the z-transform of (2.21) yields 

Sr.{^)A{z) = B{z)H{z--^). (2.22) 

Making the substitution 

m^-') = IP (2.23) 

in (2.22) and rearranging, results in 

/l(z-')S„(z)/l(j) = (2.24) 

To utilize (2.24) to find the polynomial coefficients of B{z), which are the elements 

of the sequence 6(/), we must perform spectral factorization of the sequence result- 
ing from the convolution of the three sequences a(/), a(—l), and A detailed 

explanation of this technique can be found in [Ref. 39] or [Ref. 40]. 

b. Durbin’s Method 

Durbin’s method [Ref. 41] makes use of the property by which a 
process containing zeros can be represented by an all-pole system if enough poles 
are used. The first step is to filter the sequence to be modeled through the inverse 
of the previously determined all-pole filter coefficients as shown in Figure 2.4. The 
resulting residual sequence, s(n), will nominally be an all-zero sequence. A large 
order all-pole model, Aiarffe(^), can then be fitted to s(n) to obtain the relationship 
illustrated in Figure 2.5. If the model order for Aiarffe(^) is sufficiently large, all of 
the ‘information’ in s(n) will be contained in the coefficients of Aiarffe(^)- If an all 
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x(n) 




s(n) 



Figure 2.4: Filtering process to generate the all-zero residual sequence for 
the application of Durbin’s method. 



u,(n) 





1 











s(n) 



Figure 2.5: The residual sequence approximated by a large all-pole model. 



pole model, Ij B{z), is then constructed for the sequence of coefficients in Aiarge(z)^ 
the relationship obtained is 

A:.rA^) = (2.25) 

Therefore the transfer function 1/Atarge{z) of Figure 2.5 is replaced by 



B(z) 



1 

AlaTgei^^) 



(2.26) 



which yields the desired moving average (all- zero) model, 
c. Shank’s Method 

Consider the system shown in Figure 2.6 where /i>i(n) is the impulse 
response of the previously determined all-pole portion of a pole-zero model (using 
Prony’s method for example). Shank’s method [Ref. 42] is to satisfy the relationship 
of Figure 2.6 in the least squares sense. This relationship can be described by the 
matrix equation 
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hA{Q) 


flAiQ-l) ••• 


hA[0) 




hA{Q + 1) 


hA{Q) ••• 


hA{l) 




. hA{N-l) 


hA{N-2) ••• 


hA{N-l-Q) . 





x{Q) 




esiQ) 




x[Q + 1) 


+ 


cb(Q + 1) 


(2.27) 


. x{N-l) . 




.eeiN-l) . 





H^b = x + eB. (2.28) 



Equation (2.28) can be solved in the manner of (2.6) with analogous to and x 
analagous to x^. The bk coefficients can therefore be found by using the psuedoinverse. 



x{n) 




Figure 2.6: System for Shank’s method determination of transfer function 
numerator coefficients. 

3. Iterative Prefiltering 

An iterative technique for solving the direct modeling problem of Figure 2.1 
called iterative prefiltering has been proposed in [Ref. 24]. An effective application of 
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the technique has been reported in [Ref. 43]. The presentation below follows that of 
[Ref. 40]. 

In this method, the direct modeling problem error (see Figure 2.1a), 

ed{n) = x(n) — h{n) (2.29) 



is expressed in the alternative form 

e<i(^) = x{n) — b{n) * /i>i(n) 

= x(n) * hA{n) * a{n) — b(n) * hA{n) (2.30) 



where a(n) and b(n) are the sequences of transfer function denominator and numer- 
ator coefficients, respectively, and hA{n) is the impulse response of the AR (all-pole) 
portion of the model, i.e. 

By then making the equation error iterative (superscripts represent the index of iter- 
ation), 

= x(n) * h\{n) * a‘‘*’^(n) — 6*'^^(n) * h\{n), (2.31) 



the least squares error solution for a’"*'^(n) and 6‘‘''^(n) at each iteration can be cal- 
culated using parameter estimates from the previous iteration. In matrix form (2.31) 
becomes 



^k(p) 

Xh(,P + 1 ) 



xh(N - 1) 



^i(O) h\{P) 

4(1) a;,(p+i) 

xiiN-l-P) h\{N-l) 



■ ■ • h\(P - Q) 
■■■ h\{P-Q + l) 



■■■ hU(N-l-Q)\ 



r ] 

. ei+'(W - 1) . 



J+i 

Cl 



a‘p+' 



/>‘+* 



(2.32) 
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where Xh{n) = x{n) * h\{n). Note that (2.32) can be solved in exactly the same 
manner as (2.6). 

a. Correlation Domain Iterative Prefiltering 

Many pole-zero modeling algorithm’s which were originally conceived 
based of the time domain pole-zero difference equation, (2.1), have been reformulated 
in the correlation domain. Examples of this include the correlation domain extension 
of Prony’s method resulting in the modified Yule- Walker equation methods, and an 
instrumental variable method of least squares system identification which Soderstrom 
has indicated is simply a correlation domain formulation of least squares system 
identification [Ref. 44]. In modeling trials conducted for this thesis both of these 
correlation domain methods were found to be significant improvements over their 
time domain counterparts. 

Iterative prefiltering can also be extended into the correlation domain. 
To see how this is done first note that the direct formulation of the pole-zero modeling 
problem of Figure (2.1a) may be reinterpreted in the correlation domain by employing 
the relationship of (2.20) and Figure 2.3. Figure 2.7 illustrates this new direct pole- 
zero modeling interpretation. 

Proceeding as for the time domaiin iterative prefiltering above we write 
the correlation domain error equation, 

®d,corr(^) ~ ^xx(^) ~ ®( ^) * ^(n) 

= rxx{n) * h^in) * a{n) — x{—n) * b{n) * hA{n) (2.33) 

where, as before, a{n) and b{n) are sequences of the model coefficients, hA{n) is the 
impulse response of the AR (all-pole) portion of the estimated model, and rji(n) 
is found using x{n) as the desired impulse response so that Txx(^) = x{n) * x{—n). 
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rxi(tt) 




Figure 2.7: The direct pole-zero modeling problem formulation expressed 
in the correlation domain. 

When the error is made iterative, (2.33) becomes 

* 3:(-n) * h\{n) * a‘+^(n) - x(-n) * 6’+^(n) * h\{n) (2.34) 

where the superscripts represent the index of iteration. In matrix form (2.34) becomes 

Th(P) ■■■ ri(0) x\(P) x\(P-Q) 

’•k(^’ + l) ’•■(1) x^(P+l) ■■■ xi,(P-0+l) 

_Tk{2N-l) rl(2Af-l-/>) x^(2AT-l) x^(2Af - 1 - C?) . 

r pi'X^iP) ■ 

. - 1 ) _ 

where r/,(n) = * h\{n) and Xh.{n) = x(— n) * h\{n) and which can be solved as 

before. 












(2.35) 
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4. Maximum Likelihood Techniques 

Maximum likelihood estimation of parameters is the statistical standard 
against which the performance of other estimators is measured. This estimator makes 
use of all the useful statistical information available in a given set of data [Ref. 45, p. 
73]. Difficulties arise, however, in the implementation of the maiximum likelihood es- 
timator (MLE). True maximum likelihood estimation requires exact knowledge of the 
conditional probability density function (PDF) of the observed data conditioned on 
the parameters to be estimated. This conditional PDF must then be simultaneously 
maximized with respect to all parameters being estimated. In practice, most efforts 
to employ maximum likelihood estimation meike simplifying assumptions about the 
nature of the input data to derive a useful algorithm. Such techniques are usually 
called approximate maximum likelihood methods. 

The approximate MLE chosen for this thesis is due to Akaike [Ref. 31]. 
The brief developement of this algorithm provided below follows that of Kay [Ref. 4, 
Ch. 9,10]. Additional backround on maximum likelihood estimation can be found in 
[Ref. 1,8, 32, 45] and references therein. 

Given a sequence of independent random variable observations, i(n), and 
a corresponding set of parameters to be estimated, 6k, the desired set of estimates for 
the $k's is the one for which the observed data set is the most likely. In terms of the 
conditional probability density or likelihood function, 

p(z(0), x(l), . . • , x{N - 1)1^1 , ^ 2 , • • • , Ok), 

the desired set of estimates is the one which, for a given set of x(n), is maximized. 
In the case of pole-zero modeling, an expression for the observed data sequence’s 
joint probability density function conditioned on the model peirameters, a,- and bj, is 
required. To obtain such an expression, it is generally assumed that the observed data 
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is Gaussian. If, further, it is assumed that the input to the pole zero model is white 
Gaussian noise and that the data record length is much longer than the transient 
response due to initial conditions, the conditional PDF for the observed data can be 
arrived at fairly directly. 

The joint PDF for the zero mean white Gaussian input sequence, u(n), is 
of the form 

p(u(0), u(l), • • • , u{N - 1)) = n (2-36) 

i=o \J2Tal 

where < 7 ^ is the variance. Now the density function for x(0),x(l), • • • ,x(7V — 1 ) con- 
ditioned on the can be found from (2.36) through the standard linear transfor- 
mation, 

p(x(0), x(l), • • • , x{N - 1 )) = p(u/(0), u/(l), • • • , u,{N - 1 )) \J\ , (2.37) 

where Uf{n) is the inverse filter relationship of the original pole-zero difference equa- 
tion, ( 2 . 1 ). Specifically 

1 P 1 Q 

“/(") = T “ 0 - T H - k) (2.38) 

.=0 

and is the Jacobian of the linear transformation uj. To simplify the transformation 
assume that the pole- zero filter in (2.38) has been expressed as its equivalent all-pole 
filter. 

Pap 

= XI dAP,ix{n — i). (2.39) 

i=0 

The final result of the linear transformation (2.37) will then be 

p(x(0), x(l), • • • , x{N - 1)) = n -r , exp{- ^'‘=° ^^. 40 ) 

n=o y2ircr^ 

When (2.40) is expressed in terms of the pole-zero parameters, 

oi, 02, ,ap, bo, 61 , • • • , 6q, 
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for the purpose of maximization the relationship is highly nonlinear. Note that the 
maximization of (2.40) requires the minimization, over all possible a,’s and bk's, of 



the inverse filter error, Uf{n). 

Akaike employed the Newton-Raphson iterative method for minimizing 
(2.40). This method requires the computation of Gradient and the Hessian of (2.38) 



at each iteration to generate the estimate updates 



^k+l 

















d^Q 

dadk^ 

d^Q 

dhda'^ 



d^Q 

dadh'^ 



d^Q 

dh 



(2.41) 



Using frequency domain arguments, Akaike was able to provide expressions for the 
above partial derivatives in terms of Fourier transforms which can in turn be expressed 



in terms of linear filtering operations. 

Note that there are several key assumptions made above which must be 



valid for this method of approximate MLE to apply: 



1. The data are real, Gaussian, and zero mean. 

2. The data record is large. This is to avoid end effects of assuming all data values 
outside the data record are zero when filtering the data. 

3. The poles and zeros are not close to the unit circle. This is to avoid long 
transients due to the initial conditions which are ignored. (They are assumed 
to be known and are set equal to zero.) 

At first glance these assumptions would seem to indicate that this method is inappro- 
priate to transient modeling. However, inverse filter error retains its meaning when 
considering transient waveforms; the ideal inverse filtering result for a transient signal 
is a single impulse rather than the minimum variance random sequence expected for 
a stochastic process. In fact, this is exactly the appproach taken by Jackson [Ref. 26, 
pp. 276-278] in extending Judell’s maximum likelihood method [Ref. 33] to impulse 
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response data. While this extension seems rather ad hoc, we will find that such ap- 
proximate maximum likelihood methods can be effective in modeling transient signals 
as impulse responses. A key limitation imposed by dealing with deterministic data 
is that reliance on inverse filter error excludes signals that must be modeled by non- 
minimum phase systems. In contrast, the restrictions of long data record and weak 
poles and zeros (not near the unit circle) no longer apply. Data records end effects 
and initial condition transient effects should have no impact since the assumption of 
zero valued data outside the range of data is correct if the data is chosen to end after 
most of the energy of the transient is dissipated. 

C. IMPLEMENTATION 

All modeling algorithms were implemented using the interactive language PRO- 
MATLAB from The Mathworks, Inc. on SUN workstations except for the Akaike 
MLE algorithm which was implemented in FORTRAN using a program adapted 
from [Ref. 4, Ch. 10]. The FORTRAN program was also implemented on a SUN 
workstation with a FORTRAN 77 compiler. All graphics were generated in PRO- 
MATLAB. 



23 



III. MODELING PERFORMANCE 



A, PERFORMANCE CONSIDERATIONS 

Each of the pole-zero modeling techniques presented in Chapter II is effective 
when modeling a signal that is truly the impulse response of a pole-zero system with 
no noise present and with the system order known. However, real world transient 
data rarely possess such characteristics. Real world signals of all types are notoriously 
uncooperative in fitting the signal models proposed to describe them. Reasons that 
this may be true for transient data include: 

1. Inappropriate selection of model type or modeling algorithm. 

• Linear versus non-linear models. 

• Minimum phase versus non- minimum phase rational models. 

2. The transient is time shifted because of and inappropriate selection of the data 
record starting point due to the presence of noise. 

3. The assumption of impulse system excitation is a poor approximation. 

4. Incorrect selection of model order. 

5. Noise is present in the signal. 

The test sequences used in this chapter are all obtained as the impulse response of 
linear pole-zero systems, therefore, the question of linear versus non-linear model type 
will not be at issue. For the other problems, effective transient modeling requires both 
selecting the appropriate algorithm and understanding how to use that algorithm to 
its greatest advantage. The next section describes the test sequences used in this 
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chapter. Subsequent sections address how difficulties encountered in the pole-zero 
modeling of impulse response data can arise and how they can be overcome. 

The effects of data selection, non-minimum phase systems, non-impulse excita- 
tion, and incorrect model order on modeling will initially be considered for signals 
observed with no added noise present. The effect of modeling a signal in which addi- 
tive noise is present is considered separately. We will see that situations which modify 
a linear pole-zero system often lead to another pole-zero system. This new system 
usually has the same number of poles in the same locations but with different and 
possibly additional zeros. 

B. TEST SEQUENCES 

The test impulse response sequences are generated using pole-zero models taken 
from Kay [Ref. 4]. The test sequence ARMAS uses one of Kay’s models directly while 
the test sequences ARMA4 LF, ARMAS NM, and ARMA4 CL are from Kay models 
which have been modified to enhance the illustration of certain points. The unit 
impulse response and pole-zero plot of each test sequence model is shown in Figures 
S.la-h. The model coefficients for these sequences are listed in Table S.l. 



Model 


Model Coefficients (co = = 1 ) 




Cl 


02 


<^3 


0,4 


bi 


62 


ARMAS 


-2.760 


S.809 


-2.654 


0.924 


-0.900 


0.810 


ARMA4 LF 


-S.0S5 


4.002 


-2.727 


0.778 


- 0.200 


0.040 


ARMAS NM 


-2.760 


S.809 


-2.654 


0.924 


-4.818 


25.000 


ARMA4 CL 


-2.678 


S.700 


-2.56S4 


0.917 


- 0.200 


0.040 



TABLE 3.1: Table of test sequence coefficients. 
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(a) 



(b) 



Figure 3.1: The test sequence ARMAS, (a) Impulse response plot and (b) 
pole-zero plot. 
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X model poles 




(d) 



Figure 3.1: continued The test sequence ARMA4 LF. (c) Impulse response 
plot and (d) pole-zero plot. 
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X model poles 
o model zeros 




-1.5 I 1 

z-plane 

Figure 3.1: continued The test sequence ARM A3 NM. (e) Impulse re 
sponse plot and (f) pole-zero plot. 
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z-plane 



Figure 3.1: continued The test sequence ARMA4 CL. (g) Impulse response 
plot and (h) pole-zero plot. 
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C. ESTIMATING NUMERATOR COEFFICIENTS 



1. Data Selection — Time Shifts and Initial Conditions 

Defining the range of data to be used in modeling is an important and 
usually staightforward exercise. The assumption that a signal represents the impulse 
response of a linear pole-zero system, however, implies some very specific properties 
about the initial few points of that signal. Each method of modeling the transfer 
function numerator coefficients in Chapter II reacts differently when the beginning 
points of the impulse response being modeled are degraded. Because real world 
transients do not usually exhibit the instantaneous rise time of an ideal impulse 
response and because noise is usually present, choosing the precise data range for a 
transient such as that illustrated in Figure 3.2 is often a very uncertain task. Two 
possible outcomes when the transient starting point is chosen incorrectly are: 

1. The starting point is chosen before the signal begins so that early data values 
are unrelated (and presumably of lower amplitude) to the impulse response to 
be matched (e.g. these points may consist of noise). 

2. The starting point is chosen late in which case early values of the impulse 
response are lost. 

In the first case, an adequate number of additional model zeros can account 
for the delay in the impulse response. Assunoing that the starting point for the data is 
reasonably close to the true beginning of the impulse response, any spectral features 
introduced by the unrelated early data points will probably not have enough energy to 
significantly alter the spectrum of the impulse response. Under these circumstances 
all methods can effectively find the poles. It is important to note, however, that if 
an insufficient number of zeros to account for the imposed delay is used, then some 
of the equations that are generated in Prony’s method become invalid. When these 
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Figure 3.2: An example of a laboratory generated acoustic transient. Note 
the difficulty in determining a precise starting point for the transient. 



equations are solved, the invalid equations can drastically degrade pole estimates. To 
see this we can apply Prony’s method to a system with the true orders P = 4 eind 
Q — 2. Assume the signal is delayed by inserting three zeros at the beginning of the 
data so that the original point x(0) is now the fourth data point. If we choose P = 4 
and Q = 3 in constructing the data matrix of (2.2), the resulting set of equations will 
be 



0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

x(0) 0 0 0 0 

x(l) x(0) 000 

x(2) x(l) x(0) 0 0 

x(3) x(2) x(l) x(0) 0 

x(4) x(3) x(2) x(l) x(0) 
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Note that when the lower partition is taken to find the Ck coefficients, the first two 
equations of the lower partition are invalid. Compounding the problem is that the 
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invalid equations occur in the high energy portion of a transient signal. To overcome 
this effect a numerator order of at least five is required. 

In finding the numerator coefficients of a delayed impulse response one 
of three outcomes is possible depending on the estimation technique chosen. First, 
any technique which depends directly on the initial values of the data sequence (for 
example (2.5)) will be ineffective. Second, methods which rely on the autocorrelation 
of the residual sequence will result in approximately the true zeros of the system under 
study. The original time series will not be matched directly. This case is illustrated 
in Figure 3.3. The resulting impulse response is an undelayed version of the signal. 
Finally, signal matching techniques, iterative prefiltering and Shank’s method, result 
in zeros not related to the original undelayed model but which provide the best overall 
time domain match of the delayed signal. This is shown in Figure 3.4. 

The case of choosing the data record to far to the right and thus trun- 
cating the first points of an impulse response will again have little effect on pole 
estimation. This situation corresponds the same system with (non-zero) initial con- 
ditions imposed. Since initial conditions are acounted for in the numerator, the zeros 
are significantly altered. Here the previous discussion regarding finding the under- 
lying model zeros versus obtaining a good match in the time domain still applies 
with one exception: the direct calculation of the bk coefficients from (2.5) will now 
be effective. 

2. Non-minimum Phase Modeling 

In [Ref. 46] it is demonstrated that the appropriate discrete model of 
a sampled analog waveform is often represented by a transfer function with zeros 
outside the unit circle. Many of the techniques that are currently available are pur- 
posefully structured to eliminate uch non-minimum phase models. In power spectrum 
estimation, a minimum phase system with the same frequency response magnitude 
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as a non-minimum phase system results in the same estimated spectrum. Thus for 
spectrum estimation an equivalent minimum phase system is satisfactory. In fact, 
the assumption underlying stochastic modeling techniques, namely white Gaussian 
noise input, guarantees that all transfer function combinations of minimum phase 
and maximum phase zeros are equivalent. By convention, stochastic modeling tech- 
niques always choose the minimum phase model so that the important statistic of 
inverse filter error is available. In time domain based applications, however, incor- 
rectly choosing the model phase can seriously degrade system performance. Fields 
such as seismic deconvolution, channel equilization, control, and matched filter de- 
sign, generally require identification of the correct model phase [Ref. 47, 48, 49). 
Also, we will see in Chapter IV that effective modeling of real world acoustic signals 
frequently requires non-minimum phase models. 

A number of modifications to the basic stochastic model have been in- 
troduced to allow selection of the model with the correct phase. These techniques 
generally involve changing the Gaussian nature of the input noise [Ref. 50] and often 
employ higher than second order moments or cumulants [Ref. 29, 30]. However, when 
a model’s impulse response has effectively matched a signal in the time domain, the 
resulting phase is immediately known to be the correct. Allowing for the possibility 
of non-minimum phase models places significant limitations on the modeling methods 
which may be used. Techniques which rely on inverse filtering (maximum liklihood 
methods) are not applicable since the inverse of a non-minimum phase system is unsta- 
ble. Also, techniques which use correlation information to calculate the bk coefficients 
(spectral factorization and Durbin’s method) will give poor results since correlation 
data does not preserve phase information. Figure 3.5 shows the difficulty encountered 
when Durbin’s method attempts to model a non-minimum phase system. The best 
Durbin’s method can do is produce the spectrally equivalent minimum phase version 
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of a maximum phase system since the autoregressive modeling techniques on which 
it relys can only produce zeros within the unit circle. In contrast, Figure 3.6 shows 
that Shank’s method is able to find the correct model. 

3. Numerator Modeling Summary 

Table 3.2 provides a brief summary of the modeling properties of the nu- 
merator coefficient modeling techniques considered in this thesis. Since the goal in 
Chapter IV is to perform time domain modeling of the acoustic transients being 
considered, those methods which provide the best time domain match between the 
original signal and the model impulse response are preferred. 



Method 


Equation 


Non-minimum 
phase capable? 


Data Selection 


Shifted Right 
(Delayed) 


Shifted Left 
(Truncated) 


Prony, upper 
partition 


2.5 


Yes 


Not 

Usable 


Time Series 
Match 


Spectral 

Factorization 


2.24 


No 


Underlying 

Model 


Underlying 

Model 


Durbin’s 

Method 


2.25 

2.26 


No 


Underlying 

Model 


Underlying 

Model 


Shank’s 

Method 


2.28 


Yes 


Time Series 
Match 


Time Series 
Match 


Iterative 

Prefiltering 


2.32 


Yes 


Time Series 
Match 


Time Series 
Match 


Akaike 

MLE 


2.41 


No 


Underlying 

Model 


Time Series 
Match 



TABLE 3.2: Summary of the capabilities and limitations of numerator 
modeling methods. 



D. NON-IMPULSE EXCITATION 

In any real world system the assumption of a unit impulse input is approxi- 
mate. If the duration of the excitation waveform is small relative to the period of the 
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original signal 

model impulse response 




X true poles + estimated poles 

o true zeros * estimated zeros 




Figure 3.3: The right shifted sequence ARMA3, poles modeled using 
LSMYWE and zeros modeled using Durbin’s method, (a) Time signal 
plot and (b) pole-zero plot. Durbin’s method does not account for the 
time delay but instead finds the underlying system’s true zeros. 
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original signal 

model impulse response 




X true poles + estimated poles 

o true zeros * estimated zeros 
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z-plane 



(b) 



Figure 3.4; The right shifted sequence ARMAS, poles modeled using 
LSMYWE and zeros modeled using Shank’s method, (a) Time signal 
plot and (b) pole-zero plot. The least squares method does not find the 
underlying system’s true zeros but rather achieves the best overall time 
domain match. 



36 



original signal 

model impulse response 




X true poles + estimated poles 

o true zeros * estimated zeros 




Figure 3.5: The non-minimum phase sequence ARMAS NM, poles mod- 
eled using LSMYWE and zeros modeled using Durbin’s method, (a) Time 
signal plot and (b) pole-zero plot. Durbin’s method cannot model zeros 
outside the unit circle. 
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original signal 

model impulse response 




X true poles + estimated poles 

o true zeros * estimated zeros 




z-plane 



(a) 



(b) 



Figure 3.6: The non-minimum phase sequence ARMAS NM, poles mod- 
eled using LSMYWE and zeros modeled using Shank’s method, (a) Time 
signal plot and (b) pole-zero plot. The least squares method is effective 
at modeling zeros outside the unit circle. 
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lowest oscillation frequencies present then the assumption is justified. However, it is 
reasonable to expect this criterion will often not be satisfied. Non-impulse excitations 
which may be encountered include: 

1. A long duration baseband type pulse. 

2. An uncorrelated random train of impulses. This model is often used to account 
for reverberation (echoes) in siesmic deconvolution [Ref. 51]. 

3. A frequency swept input that sweeps through the natural frequencies of a sys- 
tem. This model is usually considered in conjunction with the starting and 
stopping of rotating machinery. 

If the system input were known, the modeling problem could be formulated as a sys- 
tem identification problem. When no information about the system input is available, 
other means must be found to deal with this problem. 

1. Baseband Pulse Excitation 

The effect of modeling a transient signal from a linear system in which the 
input is a long duration, baseband-type pulse can best be understood as filtering by a 
finite impulse response (FIR, all-zero) filter as illustrated in Figure 3.7. The spectral 
properties of the original time series are windowed by the frequency response of the 
FIR filter coefficients. For this type of pulse, the effect is that of low pass filtering. 
Thus high frequencies are attenuated relative to low frequencies. If no frequency 
component exists below the FIR filter’s cutoff frequency then the original spectrum 
is altered according to the side lobe structure of the FIR filter. 

The new model that results can be viewed as a system with the original 
model poles but with new numerator polynomial coefficients that are the result of 
convolving the FIR filter coefficients with the original numerator polynomial coeffi- 
cients. This results in a higher order polynomial, hence more zeros than were present 
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^or»p(^) 



^o6s(^) 



FIR 



Figure 3.7: One way to view a linear system excited by a baseband pulse. 

in the original model will be required. Observe how the original impulse response of 
Figure 3.7a is altered to Figure 3.7b when a nine element triangular excitation pulse 
is used. The dotted lines in Figures 3.7b,c,d indicate the modeling results obtained 
when none, two, and four extra zeros, respectively, are used in the estimated pole-zero 
model. In this case four extra zeros prove sufficient to account for the input pulse. 
The corresponding model spectrum. Figure 3.7e illustrates the attenuation caused by 
the baseband excitation. The pole-zero plot in Figure 3.7f shows that non-minimum 
phase zeros were required to achieve an effective time domain match. 

2. Random Impulse Train Excitation 

If the input to a linear system is an uncorrelated random train of impulses 
then, although the time series may be significantly different from the original impulse 
response, the autocorrelation function of the signal is theoretically unaltered except 
for a scaling factor. This is because the autocorrelation of an uncorrelated impulse 
train is a scaled unit impulse. Thus modeling methods which rely on correlation 
information should be effective. However, over a finite time interval it is unlikely that 
a random sequence will be truly uncorrelated. As the impulse train becomes correlated 
the situation will be equivalent to the baseband pulse case described above. 

3. Frequency Swept Excitation 

The output of a system excited with a frequency swept signal depends on 
the rate at which the sweep occurs. A slow sweep will result in a series of transient 
events, each at a specific resonant frequency of the system. These events can each 
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be modeled separately. When the sweep rate is rapid, all natural modes will appear 
much as if the input were an impulse except that the phase relationships of the various 
components may be changed. 

E. MODEL ORDER SELECTION 

In studies of rational modeling the issue which continues to be the most con- 
founding is that of model order selection. The proposed methods which have a sound 
theoretical basis (e.g. [Ref. 52, 53, 54]) are very difficult to actually implement. 
These methods are invariably related to maximum likelihood concepts and therefore 
rely heavily on inverse filter statistics. This implies that for model order evaluation 
the inverse filter error must be calculated over all possible model orders. Then the 
model order and inverse filter error which minimize some function of the two is se- 
lected. The case of non-minimum phase systems is even more intractable since the 
inverse filter of such systems is unstable. 

A more attractive but less understood method is to initially overdetermine a 
system and then allow the modeling algorithm to indicate the correct model order. 
A method proposed by Cadzow [Ref. 39] uses singular value decomposition to aid in 
model order selection in the denominator of all-pole and pole-zero models. Kumeresan 
and Tufts [Ref. 21] have shown that when Prony’s method is applied to exponentially 
damped sinusoids reversed in time, valid poles occur outside the unit circle and excess 
poles occur inside the unit circle. In both of these methods the denominator order 
is initially overdetermined and the modeling method then provides the correct order. 
For this thesis both of these methods were applied to a number of transients. They 
were effective when applied to the noiseless impulse responses of true pole- zero systems 
but were not robust in the presence of noise or when many narrowband components 
are present. There is no similar guidance for determining the numerator model order. 
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original signal 




Figure 3.8: The sequence ARMA4 LF with (a) unit impulse excitation 
and (b) triangular pulse input modeled with correct model orders P = A 
and Q = 2 using Prony’s method and Shank’s method. 
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Figure 3.8: continued The sequence ARMA4 LF with triangular pulse 
input (c) modeled with Prony’s method and least squares identification 
with orders P = A and Q = 4 (d) modeled with Prony’s method and Shank’s 
method orders P ■= 4 and Q = 6. 
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Figure 3.8: continued The sequence ARMA4 LF with triangular input (e) 
the original impulse response and triangular output model, order (4,6), 
spectrums and (f) the corresponding pole-zero plot. Note how the base- 
band pulse has attenuated the higher frequencies. 
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Perhaps most desirable would be a modeling technique that, when overdeter- 
mined, caused excess poles and zeros to either cancel or migrate well away from the 
unit circle (all the way to the origin ideally). Tummala [Ref. 23] has demonstrated 
some success with this concept using an iterative algorithm to solve the least squares 
identification problem. Observing the degree to which the various transient modeling 
techniques of this thesis exhibit this behavior is the approach taken in the following 
comparison. This behavior was observed by modeling the ARMAS test sequence with 
different combinations of numerator and denominator order. Table 3.3 summarizes 
the results obtained. A ‘Y’ in Table 3.3 indicates a modeling technique achived an ex- 
act time domain match between the model impulse response the modeled signal. The 
most notable negative result is that both Prony’s method and LSMYWE were ineffec- 
tive when the numerator and denominator were overdetermined. Also, excess poles 
without enough zeros causes problems for correlation domain iterative prefiltering. 
Figure 3.9 is an example of the results obtained using LSMYWE when excess poles 
and zeros are present. However, when these results were used to initialize iterative 
prefiltering the excess poles and zeros were handled effectively. 

The above results and experience gained during extensive modeling trials lead to 
the following recommended strategy when modeling complex transients about which 
very little is known; 

1. Be conservative in estimating the denominator order. Signals composed of nar- 
rowband components will generally be dominated by the few components high- 
est in energy. Even if excess zeros are required to deal with a noisy signal (see 
the next section) a small number of excess zeros will usually suffice when us- 
ing LSMYWE. Pole detection techniques such as those of CadzowfRef. 5], and 
Kumeresan, and Tufts [Ref. 21] may be helpful. 
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Method 


Correctly Modeled? (Y/N) 




P=10 


P=4 


P=10 


P=6 




Q=2 


Q=10 


Q=4 


Q=10 


Prony’s 

Method 


Y 


Y 


N 


N 


LSMYWE 


Y 


Y 


N 


N 


Iterative 

Prefiltering 


Y 


Y 


Y 


Y 


Corr Domain 
Iterative Pref 


N 


Y 


Y 


Y 


Akaike 

MLE 


Y 


Y 


Y 


Y 



TABLE 3.3: The effectiveness of thesis modeling methods on ARMA3 
with varying overdetermination of model order. ‘Y’ indicates an exact 
time domain match was achieved and ‘N’ indicates a poor time domain 
match resulted. 

2. Be expansive with estimates of numerator order. (Iterative prefiltering may be 
an exception to this. See Chapter IV.) Additional zeros can help compensate 
for mistakes made in data selection and effects of non-impulse excitation. Un- 
needed zeros will usually migrate away from the unit circle. Note that when 
using Prony’s method or LSMYWE, one numerator order can be used when 
calculating the denominator coefficients and another numerator order can be 
chosen when finding the numerator coefficients. Using Shank’s method with 
Prony’s method or LSMYWE provides the most flexibility when a model is 
attempting to account for erroneous assumptions (e.g. data selection) which 
may have been made by the user. The above recomendations are primarily 
intended for Prony’s method and LSMYWE. However, using any iterative tech- 
nique requires reasonable initial estimates which are usually arrived at using 
linear estimation techniques. 
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Figure 3,9; Pole-zero plots for the sequence ARMAS modeled with overde- 
termined model orders of P=10 and Q=10. (a) Ismywe and Durbin’s 

method results in an unstable system, (b) iterative prefiltering provides 
an effective estimation of the correct poles and zeros and an excellent time 
domain match (not shown). 
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F. NOISE PERFORMANCE 



Kay [Ref. 55] has shown that for an all-pole process, additive white noise has 
the effect of introducing zeros which migrate from the origin to the model poles as 
the signal-to-noise ratio is decreased. Alternatively, if zeros are not introduced into 
the model, poles move toward the origin as the signal-to-noise ratio is reduced [Ref. 
56). In either case the overall effect is a loss of spectral resolution. This result extends 
directly to pole-zero processes in white noise. Therefore one important feature of any 
modeling technique is it’s ability to resolve spectral components in the presence of 
noise. 

1. Rational Modeling of Noisy Data 

Several authors have recently addressed the difficulties associated with 
modeling a time series in which additive white noise is present. Kay [Ref. 57] has 
noted that when the variance of the white noise can be estimated, its effect can 
be removed from the zeroth autocorrelation lag to improve the resolution of pole 
frequencies in correlation based autoregressive modeling. Alternatively, the zeroth 
autocorrelation lag may be avoided by using (2.18), sometimes called the high order 
Yule- Walker Equations, for Q > P or hy eliminating the rows containing the zeroth 
lag in certain least squares formulations [Ref. 44]. If the resulting high order Yule- 
Walker equations yield a matrix that is well conditioned, this procedure is equivalent 
to that of Kay above [Ref. 58]. Even with good correlation estimates, however, the 
matrix Ra is not guaranteed to be positive definite (i.e., invertible). Consequently, 
R^ is more likely to be a poorly conditioned matrix. This is not a problem in the 
least squares formulation, (2.19), since the very act of choosing such a formulation 
implies we expect a solution that is approximate. 

A number of authors have noted that overdetermining the number of model 
poles can have a profound effect when modeling noisy signals [Ref. 55, 56, 21, 36, 39]. 
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The explanation for this is that the extra poles are able to model the flat noise 
spectrum of the signal, preventing that portion of the signal from interfering with 
the poles needed for narrowband modeling. Finally, singular value decomposition has 
been used to aid resolution by overcoming the conditioning problems that can result 
from pole overdetermination and high order Yule- Walker equations [Ref. 5, 21, 36, 
59, 60]. 

Little work has appeared which analyzes or demonstrates the eflFectiveness 
of iterative prefiltering in the presence of additive white noise. Stoica and Soder- 
strom [Ref. 61] have shown that iterative prefiltering will perform well for very long 
stochastic data sequences in the presence of additive white noise if reasonable initial 
estimates are available to start the iterations. They also show that for arbitrary initial 
estimates iterative prefiltering is not guaranteed to converge. Tufts and Kumeresan 
[Ref. 62] have demonstrated superior performance for approximate MLE over linear 
prediction methods when modeling sinusoids in white noise. Using the approximate 
MLE method of Box and Jenkins, Cadzow found that such a method performed com- 
parably to the overdetermined high order Yule- Walker equations for the sum of two 
second order all-pole processes in white noise but with higher variance. 

The above work suggests numerous possibilities in dealing with noisy im- 
pulse response data. The effectiveness of these techniques is evaluated empirically in 
the following section by modeling noisy test sequences. 

2. Discussion and Examples 

The examples used to illustrate modeling performance are surrunarized 
in Table 3.4. All noise modeling examples were conducted using the test sequence 
ARMA4 CL (see Table 3.1). Figure 3.10 shows that Prony’s method has difficulty 
even when the SNR is as high as 30 dB although the excess poles introduced for 
Figure 3.10b dramatically increase modeling eflFectiveness. The moderate impact of 
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insuring that Q > P is illustrated for LSMYWE in Figure 3.11 in that the higher 
frequency pole moves out toward the unit circle when the modeling order is altered 
from (4,2) to (4,4). As with Prony’s method, overestimating the number of poles 
dramatically improves modeling in Figures 3.12 and 3.13. In Figure 3.14, the two 
previous noise modeling examples are again modeled with the smallest singular value 
set to zero to yield a data matrix of rank P. The technique is highly effective in 
both cases. The effectiveness of singular value decomposition when over deter mining 
the number of poles is illustrated by the examples in Figure 3.15. This is identical 
to the result described in [Ref. 21]. Discarding the excess singular values aids both 
in resolving the narrowband components present and in reducing the variance of the 
excess poles. Figures 3.16 and 3.17 show the hazards that are possible when using 
singular value decomposition. Any number of singular values discarded other than 
the those necessary to achieve a data matrix rank equal to the transfer function 
denominator order has undesirable side effects. Discarding too few singular values 
such as in the Figure 3.16 examples will resolve the desired narrowband components 
but also lead to excess poles near or even beyond the unit circle. This result will 
at best increase the likelihood of spurious narrowband components and can, in the 
worst case, lead to an unstable model. If too many singular values are discarded 
as in Figure 3.17b, the accuracy of narrowband component estimates will be badly 
degraded. Therefore, singular value decomposition, while potentially very useful in 
modeling transient signals in noise, must be used with caution when the true model 
order is unknown (which is basically always). 

The iterative prefiltering, correlation domain iterative prefiltering, and 
Akaike MLE algorithms were initialized with several of the preceeding examples and 
the results are shown in and in Figures 3.18 and 3.19. Iterative prefiltering dramati- 
cally improved the narrowband modeling of each example with which it was initialized. 
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Figure 3.18a indicates that excess poles may cause problems by being close to the unit 
circle. However, the ability of iterative prefiltering to cancel excess zeros with excess 
poles demonstrated in this chapter may mitigate this effect. Correlation domain iter- 
ative prefiltering also improved the resolution of each noise example on which it was 
tried although its convergence was slower than and the bias of its results were higher 
than standard iterative prefiltering. Also, correlation domain iterative prefiltering 
showed a alarming tendency to place excess poles near the true poles on noisy sig- 
nals. We experienced a great deal of diflSculty in trying to model noisy signals with 
the Akaike MLE algorithm. The primary problem encountered was early algorithm 
termination due to an iteration which produced a non-minimum phase zero. Akaike 
MLE was able to improve the quality of a ‘good’ Prony estimate in Figure 3.19b. 
Finally, Figure 3.20 shows a pole-zero modeling result corresponding to the example 
in Figure 3.17a. Figure 3.20a shows the noisy signal used to during modeling and 
Figure 3.20b compares the resulting model to the original noiseless impulse response. 

G. MODELING PERFORMANCE SUMMARY 

In modeling, the more that is known about the signal to be modeled, the better 
the choice of model structure and modeling algorithm that can be made. The focus 
of this thesis is the modeling of real world transient signals about which very little is 
known and whose characteristics can vary widely. The type of model is also set by the 
basic goal of this thesis, that is, a rational linear model. Under such circumstances, 
a sensible criterion in choosing a modeling algorithm is robustness. For numerator 
modeling. Shank’s method is particularly robust in the sense that it provide the 
numerator coefficients which give the minimum norm, least squares fit based on the 
previously determined denominator coeflScients. Shank’s method is insensitive to 
time shifted data records and uncertainties in model order. Durbin’s method is also 
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effective but suffers from the limitations that only minimum phase models are possible 
and the resulting models cannot account for time shifts. Since the Akaike MLE 
algorithm requires a minimum phase initial estimate, Durbin’s method can be used 
provide an such an estimate. Equation (2.5) and spectral factorization are extremely 
sensitive to degraded signals (e.g. noisy signals) and so will not be considered further. 

Denominator modeling with Prony’s method is extremely sensitive to even small 
amounts of noise. This is only slightly less true of Akaike MLE. Also, Akaike MLE 
cannot continue to iterate if a non-minimum phase model estimate is encountered. In 
contrast, LSMYWE and iterative prefiltering are quite robust in noise and iterative 
prefiltering is particularly robust to model order overdetermination. Our conclusion 
is that for modeling the real world acoustic transients of the next chapter, LSMYWE 
with Shank’s method and iterative prefiltering are likely to perform the best on the 
acoustic transients considered in the next chapter. 
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Method(order P,Q) 
(initialization) 


SNR 
in dB 


Bias Variance of Pole 

(true minus model) Estimates 


Figure 


True Poles 




.6247 ± .7509j 
.7144 ±.6711; 






Prony (4,2) 


30 


-.0598 ±.0330; 

poor 


< lO"'* 
poor 


3.10a 


Prony (6,2) 


30 


.0283 ±.0186; 
.0176 q: .0074; 


1 1 
o o 

V V 


3.10b 


LSMYWE (4,2) 


15 


.1788 ±.1872; 
.0346^.0180; 


< 10-^ 

.0012+ .0014; 


3.11a 


LSMYWE (4,4) 


15 


. 1475 ± .1035; 
.0386 q: .0045; 


< 10-^ 

.0022+ .0019; 


3.11b 


LSMYWE (6,2) 


15 


.0053 ± .0035; 
.0027 ±.00 14; 


1 1 
o o 

V V 


3.12a 


LSMYWE (8,8) 


10 


.0040 ±.0054; 
.0002 ± .0037; 


< 10-^ 
< 10"^ 


3.12b 


LSMYWE (8,8) 


5 


.0316 ±.0497; 
.0147 ± .0212; 


.0002+ .0001; 
.0008 + .0010; 


3.13a 


LSMYWE (12,12) 


5 


.0115 ±.0077; 
.0017 ± .0063; 


< lo-* 

< lo--* 


3.13b 


Prony (4,2) SVD 


30 


-.0010 ± .0040; 
-.0016 ±.0014; 


.0001 + .0002; 
< 10-^ + .003; 


3.14a 


LSMYWE (4,4) SVD 


15 


.0017^.0005; 
-.0007 ±.0005; 


A A 
o o 


3.14b 


Prony (8,2) 


20 


-.0030 ±.0387; 
.0275 ± .0188; 


.0001 + .0001; 
.0001 + .0001; 


3.15a 


Prony (8,2) SVD 


20 


-.0105 ±.0043; 
-.0059 ±.0077; 


< 10-^ 
< 10-“* 


3.15b 


LSMYWE (8,8) SVD 


5 


-.0074 ±.0100; 
-.0070 ±.0086; 


.0006 + .0008; 
.0003 + .0009; 


3.16a 


LSMYWE (8,8) SVD 


5 


-.0012 ±.0069; 
-.0062 ± .0064; 


.0004 + .0003; 
.0004+ .0008; 


3.16b 


LSMYWE (8,8) SVD 


5 


-.0037 ±.0024; 
-.0049 ±.0022; 


.0001+ .0001; 
.0002 + .0005; 


3.17a 


LSMYWE (8,8) SVD 


5 


-.0342 ±.0526; 
-.0507 ±.0730; 


.0001 + .0002; 
< 10-^ 


3.17b 


Iterative 
Prefiltering (4,4) 


15 


-.0001 ±.0001; 
.0002 ± .0001; 


A A 

O O 
1 1 


3.18a 


Iterative 
Prefiltering (8,8) 


5 


-.0001 ±.0011; 
-.0006 ±.0002; 


< 10-^ 
< 10"^ 


3.18b 


Corr Domain (8,8) 
Iter Pref 


5 


-.0054 ±.0030; 
.0002 ± .0017; 


A A 
o o 


3.19a 


Akaike MLE (6,2) 


30 


1 1 

o o 
o o 
o o 

OO 

-H -H 

o o 
o o 
o 

O CO 


A A 
o o 


3.19b 



TABLE 3.4: Examples of pole-zero modeling performance of the impulse 
response test sequence ARMA4 CL in added noise. 
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X true poles 
. estimated poles 




X true poles 
. estimated poles 




Figure 3.10: Model pole scatter plots of twenty trials illustrating the ef- 
fectiveness of overdetermining the number of model poles using Prony’s 
method. Test sequence is ARMA4 CL with 30 dB of noise, (a) Using the 
correct model order, (4,2) and (b) Using two excess poles, model order 
( 6 , 2 ). 
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X true poles 
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(a) 



(b) 



Figure 3.11; Model pole scatter plots of twenty trials illustrating the mod- 
erate benefits of zeroth lag cc^relation compensation by choosing Q = P 
when using LSMYWE. Test sequence is ARMA4 CL with 30 dB of noise, 
(a) Using the correct model order, (4,2) and (b) model order (4,4). 



55 



X true poles 
. estimated poles 




X true poles 
. estimated poles 




Figure 3.12: Model pole scatter plots of twenty trials illustrating the 
dramatic benefits of using excess poles with LSMYWE. Test sequence 
is ARMA4c. (a) Model order (6,2) with 15 dB of added noise and (b) 
model order (8,8) with lOdB of added noise. 
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X true poles 
. estimated poles 




X true poles 
. estimated poles 




Figure 3.13: Model pole scatter plots of twenty trials illustrating the 
dramatic benefits of using excess poles with LSMYWE. Test sequence 
is ARMA4 CL. (a) Model order (8,8) with 5 dB of added noise and (b) 
model order (12,12) with 5 dB of added noise. 
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X true poles 
. estimated poles 




X true poles 




Figure 3.14: Model pole scatter plots of twenty trisJs illustrating the dra- 
matic benefits of setting the smallest singular value of the data matrix to 
zero (i.e. adjust data matrix rank to P) for Prony’s method and LSMYWE. 
Test sequence is ARMA4 CL. (a) Prony’s method for model order (4,2) 
with 30 dB of added noise and (b) LSMYWE for model order (4,4) with 
15 dB of added noise. 
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X true poles 
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X true poles 
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Figure 3.15: Model pole scatter plots of twenty trials illustrating the effect 
of adjusting the data matrix rank when using excess poles for Prony’s 
method. Test sequence is ARMA4c. (a) Model order (8,2) with 20 dB 
of added noise and no rank adjustment and (b) model order (8,2) with 
20 dB of added noise with rank adjusted to P<rue = 4 using singular value 
decomposition. 
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Figure 3.16: Model pole scatter plots of twenty trials illustrating the effect 
of adjusting the data matrix rank when using excess poles for LSMYWE. 
Test sequence is ARMA4 CL. (a) Model order (8,8) with 5 dB of added 
noise with rank adjusted to Pirue + 4 = 8 and (b) model order (8,8) with 5 
dB of added noise with rank adjusted to P<rue + 2 = 6. 
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(a) 



(b) 



Figure 3.17: Model pole scatter plots of twenty trisJs illustrating the effect 
of adjusting the data matrix rank when using excess poles for LSMYWE. 
Test sequence is ARMA4 CL. (a) Model order (8,8) with 5 dB of added 
noise with rank adjusted to Ptme = 4 and (b) model order (8,8) with 5 dB 
of added noise with rank adjusted to Ptme — 1 = 3. 
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X true poles 
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(b) 



Figure 3.18: Model pole scatter plots of twenty trials illustrating the abil- 
ity of iterative prefiltering to improve the resolution of an LSMYWE es- 
timate. In both cases the iterative prefiltering algorithm was initialized 
using LSMYWE and Durbin’s method. Test sequence is ARMA4 CL. (a) 
Model order (4,4) with 15 dB of added noise and (b) model order (8,8) 
with 5 dB of added noise. 
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X true poles 
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Figure 3.19: Model pole scatter plots of twenty trials on test sequence 
ARMA4 CL illustrating (a) the ability of correlation domain prefiltering 
to improve the resolution of an LSMYWE initial estimate, order (8,8), 5 
dB added noise and (b) the ability of Akaike MLE to improve the resolution 
of a Prony’s method initial estimate, model order (6,2), 30 dB added noise. 
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original signal 

model impulse response 




Figure 3.20: Pole-zero model of the ARMA4 CL test sequence with 5 dB 
of additive white noise using LSMYWE and the smallest singular value 
discarded. The MA part was found using least squares identification. 
This example corresponds to the example in Figure (3.15a). (a) The noisy 
sequence and (b) the estimated model impulse response and the original 
noiseless sequence. Model order (8,8). 
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IV. LINEAR MODELING OF ACOUSTIC 

TRANSIENTS 



A. ACOUSTIC DATA— GENERAL 

The laboratory generated acoustic transient data available for this study were 
generated in six different ways using ordinary laboratory objects. The data records 
and their method of generation are summarized in Table 4.1. Each transient will 
be refered to by the object used or the action performed to generate that particular 
transient such as ‘book’ or ‘slam’. The six data records modeled in this section are 
shown in Figures 4.1a-f. The range of data that was actually modeled is listed in 
Table 4.1. The sampling rate for the data is unknown but is not required since there 
is no need to infer the specific characterics of the acoustic source or the acoustic 
enviroment. 



Data Name 


Generation Technique 


Indices of Modeled Data 
(begintend) 


Book 


Dropped Book 


(55:454) 


Slam 


Slammed Book 


(15:414) 


Plate 


Struck Metal Plate 


(5:704) 


Ruler 


Book Struck with Ruler 


(501:650) 


Wood 


Clapped Wooden Blocks 


(171:320) 


Wrench 


Dropped Wrench 


(101:250) 



TABLE 4.1: Summary of Acoustic Transient Data and method of genera- 
tion. 



B. ACOUSTIC TRANSIENT MODELING RESULTS 

Iterative prefiltering and correlation domain iterative prefiltering were found to 
yield the most effective time domain match (in terms of squared error) between the 
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original signal and the impulse response of the estimated model. LSMYWE with 
the smallest singular value set to zero, as in Figures 3.13a,b, was used to find the 
initializing model poles for the iterative algorithms. Removing the smallest singular 
value allowed the LSMYWE algorithm to place zeros much closer to the unit circle 
than was possible when all singular values were used. Prony’s method and the Akaike 
MLE algorithm were also used to model laboratory transients. However, their perfor- 
mance was generally poor and they will not appear in the remainder of this chapter. 
As noted in Chapter III, the problem with Prony’s method is that it cannot easily 
model highly resonant frequencies, that is, zeros close to the unit circle, in the pres- 
ence of even small amounts of noise unless a large number of excess poles are used 
and numerous singular values are discarded. This procedure is burdensome when 
initializing iterative algorithms which do not require these excess poles. The primary 
difficulty with Akaike MLE is its inability to handle zeros outside the unit circle. The 
initializing model zeros were found using Shank’s method which was by far the most 
robust algorithm for finding numerator coefficients of the methods tested in Chapter 
Three. In addition to the time domain plots of model impulse response, a pole-zero 
model spectrum and pole-zero plot is be provided for each transient so that their 
differing characteristics can be observed. All spectra were generated by squaring the 
FFT magnitude of either the model coefficients or the signal being modeled. Model 
order was chosen based on the best educated guess of the author in accordance with 
the recommendations in Chapter III, and augmented by trial and error. 

The LSMYWE model used to initialize the iterative prefiltering algorithm for 
the Slam transient is shown in Figures 4.2a and 4.3a. The corresponding iterative 
prefiltering model shown in Figures 4.2b and 4.3b. Figure 4.4a shows the iterative 
prefiltering model spectrum. Figure 4.4b illustrates a characteristic of the itera- 
tive prefiltering algorithm that was observed throughout the modeling of acoustic 
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transients. Namely, when an excess number of model parameters are used, the error 
between the model impulse response and the signal being modeled increases dramat- 
ically, mainly at the beginning of the signal. The subsequent application of Shank’s 
method is effective in reducing this initial error. This effect is shown in Figure 4.5a. 
In general, the application of Shank’s method as the last step in the modeling process 
was found to reduce the the mean squared error between the original signal and the 
model impulse response to some degree for all modeling methods. The sensitivity to 
excess parameters shown by iterative prefiltering does not affect correlation domain 
iterative prefiltering. Note that excess model zeros cause no difficulties in Figure 
4.5b. For the Book transient, two modeling trials are shown. Figures 4.6 and 4.7 
show the best time domain match obtained using iterative prefiltering and also an 
LSMYWE, Shank’s method, respectively. Although the LSMYWE model does not 
achieve as effective an impulse response match as iterative prefiltering, with SVD it 
is more sensitive to the low energy, high frequency component present in the Book 
transient at approximately The best model of Ruler is shown in Figure 4.8. The 
model spectra for Book and Ruler appear in Figure 4.9. 

The assumption that the signal being modeled is a system impulse response 
is problematic for the Plate, Wood, and Wrench signals since they do not exhibit 
the rapid decay usually associated with an impulse response. However, it is still 
possible to achieve a resonable time domain match over small segments of each signal. 
This result is illustrated in the models of in Figures 4.10, 4.11, 4.13, and 4.14. The 
model spectra for Wood and Wrench are shown in Figure 4.15. The reason that 
correlation domain iterative prefiltering was used for the Plate, Wood, and Wrench 
signals instead of standard iterative prefiltering can be illustrated by comparing the 
two techniques on a segment of the Plate signal. Although the impulse response error 
of the two models in Figures 4.10 and 4.11 are nearly identical. Figure 4.12 shows 
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that correlation domain iterative prefiltering clearly outperformed standard iterative 
prefiltering in reproducing the spectrum of the original sequence. In fact, at the more 
suitable model order of (16,16), iterative prefiltering would not converge but instead 
oscillated in a region of convergence for any model order over (12,12). Figure 4.16b 
shows the model impulse response obtained when the large segment of Plate shown 
in Figure 4.16a is modeled using iterative prefiltering. Although the time domain and 
spectral properties of the model relative to the original signal are considerable poorer 
than those obtained for a short segment, the model does clearly share many of the 
features of the original signal. 

C. ACOUSTIC SIGNAL MODELING SUMMARY 

The modeling results obtained in the previous section of this Chapter indicate 
the possible utility of pole-zero modeling algorithms with regard to modeling tran- 
sient signals. Signals with decaying narrowband components (e.g. Slam, Book, and 
Ruler) and signals with substained narrowband components (e.g. Plate, Wood, and 
Wrench) can be modeled as the impulse response of a rational linear system. Robust 
modeling algorithms are available which can effectively deal with the many uncertain- 
ties associated with real world signals. Although the goal of achieving an exact time 
domain match between the original signal and the pole-zero model impulse response 
was not realized for any of the acoustic signals in this chapter, in all cases the degree 
of match obtained clearly indicates that many signal characteristics are described by 
the model. 



68 





Figure 4.1: Laboratory generated acoustic transient data: (a) Slam and 
(b) Book. 
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Figure 4.1; continued Laboratory generated acoustic transient data: (c) 
Ruler and (d) Plate. 
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Figure 4.1: continued Laboratory generated acoustic transient data: (e) 
Wood and (f) Wrench. 
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Figure 4.2; Modeling the transient Slam — model impulse response versus 
the original signal. Model order (6,8). (a) LSMYWE and Shank’s method 
and (b) iterative prefiltering initialized with (a). 
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Figure 4.3: Modeling the transient Slam — model pole-zero plots. Model 
order (6,8). (a) LSMYWE with SVD and Shank’s method and (b) iterative 
prefiltering initialized with (a). 
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Figure 4.4; Modeling the transient Slam — model spectrum and and over- 
parametrized iterative prefiltering model impulse response versus the orig- 
inal signal, (a) Model spectrum for iterative prefiltering with model order 
(6,8) and (b) iterative prefiltering for model order(6,12). Note how excess 
model parameters cause error at the beginning of the iterative prefiltering 
model. 
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Figure 4.5: Modeling the transient Slam — model impulse response versus 
the original signal, overparameterized modeling effects, (a) The applica- 
tion of Shanks method to the iterative prefiltering model of order (6,12) 
reduces the error at the beginning of the model impulse response and 
(b) correlation domain iterative prefiltering for model order (6,12) can 
effectively use the excess model parameters. 
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Figure 4.6: Modeling the transient Book — model impulse response versus 
the original signal.. Model order (6,6). (a) LSMYWE with SVD and 

Shank’s method and (b) iterative prefiltering initialized with (a). Note 
the sensitivity of LSMYWE with SVD to the high frequency components 
present. 
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Figure 4.7: Modeling the transient Book — model pole-zero plots. Model 
order (6,6). (a) LSMYWE with SVD and Shank’s method and (b) iterative 
prefiltering initialized with (a). Note the sensitivity of LSMYWE with 
SVD to the high frequency pole pair. 
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Figure 4.8: Modeling the transient Ruler — Model order (6,16). (a) Iter- 
ative prefiltering model impulse response versus the original signal and 
(b) the corresponding model pole-zero plot. Note that the two beating 
sinusoids have been effectively modeled. 
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Figure 4.9: Modeling the transient Book and Ruler — model spectrums. 
(a) The iterative prefiltering model of order (6,6) for the transient Book 
and (b) the iterative prefiltering model of order (6, 12) for the transient 
Ruler. 
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Figure 4.10: Modeling the transient Plate — Model order (12,12). (a) It- 
erative prefiltering model impulse response versus the original signal and 
(b) the corresponding model pole-zero plot. 
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Figure 4.11: Modeling the transient Plate — Model order (16,16). (a) Cor- 
relation domain iterative prefiltering model impulse response versus the 
original signal and (b) the corresponding model pole-zero plot. 
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Figure 4.12: Modeling the transient Plate — model spectrums versus the 
original signal spectrums. (a) The iterative prefiltering model of order 
(12,12) and (b) the correlation domain iterative prefiltering model of order 
(16,16) 
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Figure 4.13: Modeling the transient Wood — ^Model order (16,16). (a) 

Correlation domain iterative prefiltering model impulse response versus 
the original signal and (b) the corresponding model pole-zero plot. 
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Figure 4.14: Modeling the transient Wrench — Model order (8,8). (a) Cor- 
relation domain iterative prefiltering model impulse response versus the 
original signal and (b) the corresponding model pole-zero plot. 
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Figure 4.15: Modeling the transient Wood and Wrench — model spectrums. 
(a) The iterative prefiltering model of order (16,16) for the transient Wood 
and (b) the iterative prefiltering model of order(8,8) for the transient 
Wrench. 
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(b) 



r4?n i 9 «;n\ '• the transient Plate over the large segment 

( 0^250) using Iterative prefiltering, model order (16,16). (a) The orig- 
inal Plate segment and (b) the model impulse response. ® 
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V. CONCLUSIONS 



A. PERFORMANCE COMPARISON SUMMARY 

The modeling of signals as the impulse response of a linear pole-zero system is 
an important tool in digital signal processing. One natural area for applying such 
an approach is in the modeling of transient, impulse response-like waveforms. The 
specific approach taken in this thesis was to determine which pole- zero modeling al- 
gorithms are most suited to modeling complex, real world transient waveforms. The 
modeling criterion emphasized is to obtain the best (least squares error) time domain 
match between the model impulse response and the original signal. This criterion 
was adopted because it provides a degree of signal characterization that is a step be- 
yond normal power spectrum estimation. Indeed, the strength of pole-zero models is 
their ability to describe not only the resonances present (model poles), but also how 
these resonances are related (model zeros). Because of the widely varying character- 
istics anticipated for real world signals, a key evaluation criterion for any modeling 
algorithm is robustness. In particular, algorithms must be effective in the presence of 
signal degrading effects like noise and model degrading effects such as unknown model 
order. The performance of several selected algorithms were compared for known im- 
pulse response test sequences in Chapter III. The modeling experience gained in these 
experiments was then applied to modeling laboratory generated acoustic transient 
datain in Chapter IV as a test of ‘real world’ effectiveness. 

Four basic algorithms were chosen for comparison: Prony’s method, the least 
squares modified Yule- Walker equations (LSMYWE), iterative prefiltering, and the 
Akaike maximum likelihood estimator (MLE). An algorithm which is an extension of 
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iterative prefiltering into the correlation domain was also presented. For those meth- 
ods in which the model poles and zeros are determined separately (Prony’s method 
and LSMYWE), four methods for determining model zeros (i.e., transfer function nu- 
merator coefficients) were considered: the upper partition of Prony’s method, spectral 
factorization, Durbin’s method, and Shank’s method. The major conclusions of the 
algorithm performance comparisons conducted in Chapter III and Chapter IV are as 
follows: 

1. Algorithms that are unable to model zeros outside the unit circle (Durbin’s 
method, Akaike MLE) have limited versatility when modeling arbitrary tran- 
sient waveforms. All the acoustic transients in Chapter IV required a non- 
minimum phase model to obtain the best time domain match. 

2. The most robust and effective method for finding zeros that gives the best least 
squares time domain match was found to be Shank’s method. In fact, applying 
Shank’s method as a last step improved the final model of all algorithms to 
some degree. The upper partition of Prony’s method and spectral factorization 
are not very useful because of their extreme sensitivity to noisy or time shifted 
signals. 

3. Prony’s method and Akaike MLE have difficulty modeling signals in which 
additive noise is present. Even small amounts of additive noise causes a dramatic 
loss of pole frequency resolution for Prony’s method. The use of excess poles 
and singular value decomposition were found to be effective in overcoming these 
effects but these methods depend on a knowledge of correct model order. Also, 
unlike spectrum estimation, excess poles must be retained for time domain 
matching. This is in direct contrast to the parsimonious use of model parameters 
normally provided by a pole-zero model. The difficulties encountered with the 
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Akaike MLE were computational in nature; for noisy signals the Akaike MLE 
algorithm usually terminated prematurely when either a non- minimum phase 
model was encountered during an iteration or when a numerical overflow was 
induced by an unstable model estimate. 

4. LSMYWE with singular value decomposition and iterative prefiltering (includ- 
ing the correlation domain version) were found to be the most effective algo- 
rithms for modeling a signal when additive noise is present. If the true model 
order of a system is unknown, it is best to discard only one singular value. 
Almost all the resolution gain occurs with the first singular value. Discarding 
additional singular values is intended to reduce the variance of excess poles but 
it will cause poles to be biased if all model poles are necessary for signal model- 
ing. Both of these methods demonstrated the consistent ability to model poles 
very close to the unit circle. This capability was essential when modeling the 
acoustic transients used in this thesis. 

Combining the above observations leads to our recommended strategy for mod- 
eling an arbitrary transient signal. First, use LSMYWE with one singular value re- 
moved and Shank’s method to find an initial model estimate. Next, use this estimate 
to initialize either iterative prefiltering or correlation domain iterative prefiltering. Fi- 
nally, if desired, apply Shank’s method to optimize the time domain fit of the model 
zeros. 

B. RECOMMENDATIONS FOR FUTURE STUDY 

The results of Chapter IV indicate that there are pole-zero modeling algorithms 
available that are sufficiently robust to be useful for modeling many complex, real 
world transients. A number of issues regarding the pole-zero modeling of transient 
signals and the application of such models require further study. These issues include: 
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1 . The degree to which the noise performance of correlation domain iterative pre- 
filtering may be improved by introducing noise compensation of the autocorre- 
lation sequence requires study. 

2. A study of the effectiveness of model order determination techniques when ap- 
plied to transient modeling would facilitate more effective use of transient mod- 
eling techniques. 

3. The effectiveness of iterative prefiltering and correlation domain prefiltering as 
a spectral estimation technique should be explored further. 

4. The ability of pole-zero models to describe the time domain characteristics of 
a signal could aid in the detection of signals. Such an application needs to be 
persued. 

5. A study of the structural relationship of pole- zero models to the specific systems 
that generate transients may increase the usefulness of such models. 



90 



REFERENCES 



1. George E. P. Box and Gwilyn M. Jenkins, Time series analysis, forecasting and 
control, rev. ed., Holden Day, San Francisco, 1976. 

2. John M. Gottman, Time-series analysis, a comprehensive introduction for social 
scientists, Cambridge University Press, 1981. 

3. Mike West and Jeff Harrison, Bayesian forecasting and dynamic models. Springer 
Verlag, New York, 1989. 

4. Steven M. Kay, Modern spectral estimation: theory and applications, Prentice 
Hall Inc., Englewood Cliffs, New Jersey, 1988. 

5. James A. Cadzow, “Spectral estimation; an overdetermined rational model equa- 
tion approach,” Proc. IEEE, Vol. 70, No. 9, September 1982. 

6. S.B. Kesler, ed., ’’Modern spectrum analysis II,” selected reprint series of IEEE 
Acoust. Speech and Sig. Proc. Society, IEEE Press, 1986. 

7. N.K. Sinha and B. Kuszta, ’’Modeling and identification of dynamic systems,” 
Van Nostrand Reinhold Company, New York, 1983. 

8. Pieter Eykoff, “System identification,” John Wiley and Sons, New York, 1974. 

9. Lennart Ljung, System identification: theory for the user, Prentice-Hall Inc., 
Englewood Cliffs, New Jersey, 1987. 

10. Torsten Soderstrum and Petre Stoica, System identification, Prentice Hall, Lon- 
don, 1989. 

11. Lawrence R. Rabiner and Ronald W. Schafer, Digital signal processing of speech 
signals, Prentice-Hall Inc., Englewood Cliffs, New Jersey, 1978. 

12. R. Genesio and M. Milanese, “A note on the derivation and use of reduced-order 
models,” IEEE Trans. Automat. Contr., Vol. AC-21, pp. 118-122, February 1976. 

13. Cyril M. Harris, ed.. Shock and vibration handbook, McGraw-Hill Book Company, 
New York, 1988. 

14. D. Brown, G. Carbon, and K. Ramsey, "Survey of excitation techniques applicable 
to the testing of automotive structures,” Soc. Automotive Engineers, SAE Paper 
77029, 1977. 

15. C.E. Baum, ’’Emerging technology for transient and broad-band analysis and 
synthesis of antenna scatterers,” Proc. IEEE, Vol. 64, pp. 1598-1616, 1976. 

16. Dale G. Stone, “Wavelet estimation,” Proceed. IEEE, Vol. bf 72, No. 10, pp. 
1394-1402, October 1984. 

17. Micheal A. Van Blaricum and Raj Mittra, ’’Problems and solutions associated 
with Prony’s method for processing transient data,” IEEE Trans, on Ant. and 
Prop., Vol. AP-28, No. 1, pp. 174-182, January 1978. 

18. T.L. Henderson, ’’Geometric methods for determining system poles from transient 
response,” IEEE Trans, on Acoust., Speech, and Sig. Proc., Vol. ASSP-29, No. 5, 
pp. 982-988, October 1881. 



91 



19. Sophocles J. Orfanides, "Pole retrieval by eigenvector methods,” IEEE Inti. Conf. 
Acoust., Speech, and Sig. Proc., pp. 1505-1508, 1987. 

20. D.L. Brown, et.al., “Parameter estimation techniques for modal analysis,” Soc. 
Automotive Engineers, SAE Paper 790221, 1979. 

21. Ramdas Kumaresan and Donald W. Tufts, ’’Estimating the parameters of ex- 
ponentially damped sinusoids and pole-zero modeling in noise,” IEEE Trans. 
Acoust. Speech Sig. Proc., Vol. ASSP-30, pp. 833-840, December 1982. 

22. James H. McClellan, “Parametric signal modeling,” in Advanced topics in signal 
processing, Jae S. Lim and Alan V. Oppenheim, eds., pp. 1-57, Prentice Hall, 
Englewood Cliffs, New Jersey, 1988. 

23. Murali Tummula, “Iterative method for ARMA least squares identification,” sub- 
mitted for publication. 

24. Steiglitz, K. and McBride, L.E., “A technique for the identification of linear 
systems,” IEEE Trans. Automat. Contr., Vol. AC-10, pp. 461-464, October 1965. 

25. Alan G. Evans and Robert Fischl, “Optimal least squares time-domain synthesis 
of recursive digital filters,” IEEE Trans. Audio Electoacoust., Vol. AU-21, No. 1, 
pp. 61-65, February 1973. 

26. Jackson, L.B., “Digital Filters and Signal Processing,” Second Edition, Kluwer 
Academic Publishers, Boston, 1989. 

27. Ira S. Konvalinka and Miroslav R. Matausek, “Simultaneous estimation of poles 
and zeros in speech analysis and ITIF — itrative inverse filtering algorithm,” IEEE 
Trans. Acoust., Speech and Signal Process., Vol. ASSP-27, No. 5, pp. 485-492, 
October 1979. 

28. M. Isabel Ribeiro and Jose M.F. Moura, “Dual algoritm for ARMA spectrum 
estimation,” IEEE Conf. Acoust., Speech and Signal Process., pp. 2081-2084, 
1987. 

29. Jitendra K. Tugnait, “Identification of linear stochastic systems via second- and 
fourth-order cumulant matching,” IEEE Trans. Inform. Theory, Vol. IT-33, No. 
3, pp. 393-407, May 1987. 

30. Georgios B. Giannakis and Jerry M Mendel, “Stochastic realization of non- 
minimum phase systems,” Proc. American Control Conf. (ACC), pp. 1254-1259, 
Seattle, WA, June 1986. 

31. H. Akaike, “Maximum likelihood identification of gaussian autoregressive moving- 
average models,” Biometrica, Vol. 60, pp. 255-265, August 1973. 

32. Craig F. Ansley, “An algorithm for the exact likelihood of a mixed autoregressive- 
moving average process,” Biometrika, Vol. 66, No. 1, pp. 59-65, 1979. 

33. N.H. Judell, “Maximum likelihood parameter estimation for signals with rational 
spectra,” M.S. Thesis, Univ. of Rhode Island, 1983. 

34. Leland B. Jackson, Jianguo Huang, and Kevin P. Richards, “AR, ARMA, and 
AR-in-noise modeling by fitting windowed correlation data,” IEEE Conf. Acoust, 
Speech and Signal Process., pp. 2039-2042, 1987. 

35. R. Moses et. al., “An eflacient linear method of ARMA spectral estimation,” IEEE 
Conf. Acoust., Speech and Signal Process., pp. 2077-2080, 1987. 



92 



36. Donald W. Tufts and Ramdas Kumaresan, “Estimation of frequencies of multi- 
ple sinusoids; making linear prediction perform like maximum likelihood,” Proc. 
IEEE, Vol. 70, No. 9, pp. 975-989, September 1982. 

37. Gerard Alengrin and Josiane Zerubia, “A method to estimate the parameters 
of an ARMA model,” IEEE Trans. Automat. Contr., Vol. AC-32, No. 12, pp. 
1113-1115, December 1987. 

38. Haykin, S., “Modern Filters,” Macmillan Publishing Company, New York, 1989. 

39. James A. Cadzow, “Spectral analysis,” in D. F. Elliot, ed.. Handbook of digital 
signal processing, Academic Press, Inc., 1987. 

40. Therrien, C.W., “Discrete random signals and statistical signal processing,” Un- 
published Manuscript, 2 April 1990. 

41. J. Durbin, “Efficient estimation of parameters in moving- average models,” 
Biometrika, Vol. 46, pp. 306-316, 1959. 

42. John L. Shanks, “Recursion filters for digital processing,” Geophysics, Vol. 32, 
No. 1, pp. 33-51, February 1967. 

43. Steiglitz, K., “On the simultaneous estimation of poles and zeros in speech anal- 
ysis,” IEEE Trans. Acoust. Speech Signal Proc., Vol. ASSP-25, pp. 229-234, June 
1977. 

44. Torsten Soderstrom, “Identification of stochastic linear systems in the presence 
of input noise,” Automatica, Vol. 17, No. 5, pp. 713-725, 1981. 

45. Sir Ronald A. Fisher, ''''Statistical methods and scientific inference,” Third Ed., 
Hafner Press, A Division of Macmillan Publishing Co., Inc., New York, 1973. 

46. K.J. Astrom, P. Hagander, and J Sternby, “Zeros of sampled systems,” Automat- 
ica, Vol. 20, pp. 31-38, 1984. 

47. J.M. Mendel, '^Optimal Siesmic Deconvolution,” Academic Press, New York, 
1983. 

48. Albert Benveniste, Maurice Goursat, and Gabriel Ruget, “Robust identification 
of a nonminimum phase system: blind adjustment of a linear equalizer in data 
communications,” IEEE Trans. Automat. Contr., Vol. AC-25, No. 3, pp. 385-398, 
June 1980. 

49. J.S. Freudenburg and D.P. Looze, “Sensitivity reduction, nonminimum phase 
zeros, and design tradeoffs in single loop feedback systems,” IEEE Inti. Conf. 
Dec. and Contr., pp. 625-630, December 1983. 

50. John J. Kormylo and Jerry M Mendel, “Identifiability of nonminimum phase 
linear stochastic systems,” IEEE Trans. Automat. Control, Vol. AC-28, No. 12, 
pp. 1081-1089, December 1983. 

51. Vijay K. Arya and H.D. Holden, “Deconvolution of seismic data — an overview,” 
IEEE Trans. Geosci. Electron., Vol. GE-16, No. 2, pp. 95-98, April 1978. 

52. Hirotugo Akaike, “A new look at statistical model identification,” IEEE Trans. 
Automat. Contr., Vol. AC-19, No. 6, pp. 716-723, December 1974. 

53. E.J. Hannan, “The estimation of the order of an ARMA process,” Annals of Stat., 
Vol. 8, No. 5, pp. 1071-1081, 1980. 



93 



54. Rangaisaoni Kashyap, “Optimal choice of AR and MA parts in autoregressive 
moving average models,” IEEE Trans. Patt. Analy. and Mach. Intelligence.! Vol. 
PAMI-4, No. 2, pp. 99-104, March 1982. 

55. Steven M.Kay, “The effects of noise on the autoregressive spectral estimator,” 
IEEE Trans, on Acoust, Speech, and Signal Proc., Vol. ASSP-27, No. 5, pp. 478- 
485, October 1979. 

56. Leland B. Jackson, et. al., “Frequency estimation by linear prediction,” Proc. 
Inti. Conf. Acoust., Speech, Signal Proc., pp. 352-356, 1978. 

57. Steven M. Kay, “Noise compensation for autoregressive spectral estimates,” IEEE 
Trans. Acoust., Speech, Signal Process., Vol. ASSP-28, No. 3, pp. 292-303, June 
1980. 

58. Luis Vergara-Dominguez, “New insights into the high-order Yule- Walker equa- 
tions,” IEEE Trans. Acoust., Speech, Signal Process., Vol. ASSP-38, No. 9, pp. 
1649-1651, September 1990. 

59. James A. Cadzow and Behshad Baseghi, “Data adaptive ARMA modeling of time 
series,” Proc. Inti. Conf. Acoust., Speech, Signal Process., pp. 256-261, 1982. 

60. Srinivasan, T., Swanson, D.C., Symons, F.W., “ARMA model order/data length 
tradeoff for specified frequency resolution,” Proc. ICASSP, 1984, Paper 38.2. 

61. Petre Stoica and Torsten Soderstrom, “The Steiglitz-McBride identification algo- 
rithm revisited — convergence analysis and accuracy aspects,” IEEE Trans. Au- 
tomat. Contr., Vol. AC-26, No. 3, pp. 712-717, June 1981. 

62. Donald W. Tufts and R. Kumaresan, “Improved spectral resolution,” Proc. IEEE, 
Vol. 68, No. 3, pp. 419-420, March 1980. 



94 



INITIAL DISTRIBUTION LIST 



No. of Copies 

1. Defense Technical Information Center 2 

Cameron Station 

Alexandria, Virginia 22304-6145 

2. Library, Code 52 2 

Naval Postgraduate School 

Monterey, California 93943-5002 

3. Chairman, Code EC 1 

Department of Electrical and 

Computer Engineering 
Naval Postgraduate School 
Monterey, California 93943-5000 

4. Superintendent, Naval Postgraduate School 7 

Atten: Professor Charles Therrien, Code EC/Th 

Naval Postgraduate School 
Monterey, California 93943-5000 

5. Superintendent, Naval Postgraduate School 1 

Atten: Professor James H. Miller, Code EC/Mr 

Naval Postgraduate School 
Monterey, California 93943-5000 

5. LT Gary L. May 1 

5545 N. 83rd Ave. 

Glendale, Arizona 85303 



95 




^9 A '^C:i 







4 











